require(list)

#### data preparation ####
FtF.data <- read.csv("FtF_data.csv")
online.data <- read.csv("online_data.csv")

## reversed-coded DQ outcome for combinedListDirect()
FtF.data$DQ.rev <- 1 - FtF.data$DQ
online.data$DQ.rev <- 1 - online.data$DQ

## reversed-coded list outcome for combinedListDirect()
FtF.data$list.rev <- ifelse(FtF.data$treatment == 1, 
                            4 - FtF.data$list, 3 - FtF.data$list)
online.data$list.rev <- ifelse(online.data$treatment == 1, 
                               4 - online.data$list, 3 - online.data$list)

#### test for the no design effect assumption (Appendix A.1) ####
## face-to-face
no.design.effect.test.FtF <- ict.test(y = FtF.data$list, 
                                      treat = FtF.data$treatment, J = 3)
round(no.design.effect.test.FtF$pi.table, 3)
no.design.effect.test.FtF$p

## online
no.design.effect.test.online <- ict.test(y = online.data$list, 
                                         treat = online.data$treatment, J = 3)
round(no.design.effect.test.online$pi.table, 3)
no.design.effect.test.online$p

#### discussion of the no liars assumption (Appendix A.2) ####
round(prop.table(table(FtF.data$treatment, FtF.data$list), margin = 1), 2)
round(prop.table(table(online.data$treatment, online.data$list), margin = 1), 2)

cairo_pdf("Figure_A1.pdf", width = 6, height = 3, pointsize = 8)
layout(matrix(1:2, 1, 2, byrow = TRUE))
par(mar = c(3, 4, 3, 1), lwd = 0.5)
hist(FtF.data$list[FtF.data$treatment == 1], breaks = 0:5, 
     freq = FALSE, right = FALSE, col = "white", 
     ylim = c(0, 0.6), axes = FALSE, 
     main = "Face-to-Face Survey", xlab = "", ylab = "Proportion")
hist(FtF.data$list[FtF.data$treatment == 0], breaks = 0:5, 
     freq = FALSE, right = FALSE, 
     col = "#A0A0A050", border = "#A0A0A0", density = 50, add = TRUE)
legend("topright", legend = c("Control", "Treatment"), 
       fill = c("#A0A0A050", "white"), border = c("#A0A0A0", "black"), 
       bty = "n", density = c(50, NA))
segments(0, 0, 4, 0)
axis(2, lwd = 0.5)
mtext(0:4, side = 1, at = 0:4 + 0.5, line = 0)
mtext("Count", side = 1, line = 1.5)
hist(online.data$list[online.data$treatment == 1], breaks = 0:5, 
     freq = FALSE, right = FALSE, col = "white", 
     ylim = c(0, 0.6), axes = FALSE, 
     main = "Online Survey", xlab = "", ylab = "Proportion")
hist(online.data$list[online.data$treatment == 0], breaks = 0:5, 
     freq = FALSE, right = FALSE, 
     col = "#A0A0A050", border = "#A0A0A0", density = 50, add = TRUE)
legend("topright", legend = c("Control", "Treatment"), 
       fill = c("#A0A0A050", "white"), border = c("#A0A0A0", "black"), 
       bty = "n", density = c(50, NA))
segments(0, 0, 4, 0)
axis(2, lwd = 0.5)
mtext(0:4, side = 1, at = 0:4 + 0.5, line = 0)
mtext("Count", side = 1, line = 1.5)
dev.off()

#### additional combined tests of several assumptions (Appendix A.3) ####
## face-to-face
Aronow.et.al.FtF <- combinedListDirect(list.rev ~ treatment, data = FtF.data, 
                                       treat = "treatment", direct = "DQ.rev")

# first test
round(1 - Aronow.et.al.FtF$placebo.I$estimate, 3)
round(Aronow.et.al.FtF$placebo.I$p, 3)

# second test
round(Aronow.et.al.FtF$placebo.II$estimate, 3)
round(Aronow.et.al.FtF$placebo.II$p, 3)

## online
Aronow.et.al.online <- combinedListDirect(list.rev ~ treatment, data = online.data, 
                                          treat = "treatment", direct = "DQ.rev")

# first test
round(1 - Aronow.et.al.online$placebo.I$estimate, 3)
round(Aronow.et.al.online$placebo.I$p, 3)

# second test
round(Aronow.et.al.online$placebo.II$estimate, 3)
round(Aronow.et.al.online$placebo.II$p, 3)